# Marijuana experiment

# Load relevant libraries
library(data.table)
library(tidyverse)
library(broom)
library(Cairo)
library(ggthemes)
library(knitr)
library(kableExtra)
library(psych)
library(xtable)

red_mit <- '#A31F34'
red_light <- '#A9606C'
blue_mit <- '#315485'

############################
# Data Processing
############################

w1 <- read_csv("Data/MediaSSI_Dec2017_w1_recoded.csv") %>%
  filter(consent == 1)
w2 <- read_csv("Data/MediaSSI_Dec2017_w2_recoded.csv")
w3 <- read_csv("Data/MediaSSI_Dec2017_w3_recoded.csv")

keep_best <- \(df, prog) df %>%
  filter(!is.na(.data[[prog]])) %>%
  arrange(PID, desc(.data[[prog]])) %>%
  group_by(PID) %>% slice_head(n = 1) %>% ungroup()

w1u <- keep_best(w1, "Progress")
w2u <- keep_best(w2, "Progress")
w3u <- keep_best(w3, "Progress")

w123 <- w1u %>%
  left_join(w2u, by = "PID", suffix = c("", "_w2")) %>%
  left_join(w3u, by = "PID", suffix = c("", "_w3"))



# Define the outcomes
outcomes <- c("mar_tradeoff", "mar_econ", "mar_costmore", "mar_fewserious",
              "mar_wrong", "mar_violence", "mar_legmed", "mar_serious",
              "mar_legrec", "danger_mar")

# Calculate Cronbach's alpha
index_fa <- psych::alpha(w123[, outcomes], check.keys = TRUE)
alpha <- index_fa$total["raw_alpha"]
print(paste("Cronbach's alpha:", round(alpha, 2)))

w123 <- w123 %>%
  mutate(
    mar_index = rowMeans(.[outcomes], na.rm = TRUE),
    mar_index_w2 = rowMeans(.[paste0(outcomes, "_w2")], na.rm = TRUE),
    mar_index_w3 = rowMeans(.[paste0(outcomes, "_w3")], na.rm = TRUE)
  )

calculate_pca <- function(data, vars) {
  result <- rep(NA, nrow(data))
  complete_cases <- complete.cases(data[vars])
  if(sum(complete_cases) > 0) {
    complete_data <- data[complete_cases, vars]
    pca <- prcomp(complete_data, center = TRUE, scale. = TRUE)
    pca_scores <- pca$x[,1]
    
    # Check correlation with the additive index
    additive_index <- rowMeans(complete_data, na.rm = TRUE)
    if(cor(pca_scores, additive_index) < 0) {
      pca_scores <- -pca_scores  # Flip the sign if negatively correlated
    }
    
    pca_scores <- scale(pca_scores)  # Standardize
    result[complete_cases] <- pca_scores
  }
  return(result)
}

w123$pcw1 <- calculate_pca(w123, outcomes)
w123$pcw2 <- calculate_pca(w123, paste0(outcomes, "_w2"))
w123$pcw3 <- calculate_pca(w123, paste0(outcomes, "_w3"))

###############################################
# Figure 3 - Persuasive effects in full sample
###############################################

extract_coefs <- function(model) {
  coefs <- coef(model)
  se <- sqrt(diag(vcov(model)))
  data.frame(
    estimate = coefs[-1],
    se = se[-1],
    treatment = c("Fox", "MSNBC"),
    conf.low = coefs[-1] - qnorm(0.975) * se[-1],
    conf.high = coefs[-1] + qnorm(0.975) * se[-1],
    conf.low.90 = coefs[-1] - qnorm(0.95) * se[-1],
    conf.high.90 = coefs[-1] + qnorm(0.95) * se[-1]
  )
}

# First Principal Component Plot
pc_w1 <- lm(pcw1 ~ fox + msnbc, data = filter(w123, forcedchoice == 1 & forcedchoice_w2 == 0))
pc_w2 <- lm(pcw2 ~ fox + msnbc, data = filter(w123, forcedchoice == 1 & forcedchoice_w2 == 0))
pc_w3 <- lm(pcw3 ~ fox + msnbc, data = filter(w123, forcedchoice == 1 & forcedchoice_w2 == 0))

fig1_data_pc <- bind_rows(
  extract_coefs(pc_w1) %>% mutate(wave = 1),
  extract_coefs(pc_w2) %>% mutate(wave = 2),
  extract_coefs(pc_w3) %>% mutate(wave = 3)
) %>% mutate(outcome = "First Principal Component")

labs_pc <- data.frame(
  wave=c(1.5, 1.5), 
  estimate=c(0.25, -0.2), 
  label=c("Fox rather\nthan Entertainment", "MSNBC rather\nthan Entertainment"),
  treatment=c("Fox", "MSNBC")
)

pdf('Output/fig3_pc.pdf', width=6, height=4)
print(ggplot(fig1_data_pc, aes(y=estimate, x=wave, color=treatment)) +
  geom_hline(yintercept = 0, lty=2) + 
  geom_point(position=position_dodge(width=0.25))+ 
  geom_line(aes(group=treatment), position=position_dodge(width=0.25)) + 
  geom_errorbar(aes(ymin=conf.low, ymax=conf.high, width=0), size=0.5, position=position_dodge(width=0.25)) + 
  geom_errorbar(aes(ymin=conf.low.90, ymax=conf.high.90, width=0), size=1, position=position_dodge(width=0.25)) + 
  geom_text(data=labs_pc, aes(x=wave, y=estimate, label=label, color=treatment)) +
  theme_bw() + 
  scale_x_continuous('Survey Wave', breaks=c(1,2,3), limits=c(0.75,3.25)) + 
  scale_y_continuous('Effect of Treatment on\nFirst Principal Component',
                     breaks=seq(-0.5,0.5,0.25),
                     labels=c("\n\n-0.5\n(More\nliberal)","-0.25",  "0.0","0.25", "(More\nconservative)\n0.5\n\n"),
                     limits=c(-0.55,0.55)) + 
  scale_color_manual("Treatment", breaks = c("Fox","MSNBC"), values=c(red_mit,blue_mit)) +
  theme(legend.position="none", strip.background = element_rect(fill="lightgrey"), axis.title.y = element_text(margin = margin(r = -20))))
dev.off()


## Figure 3 notes

# Analysis sample
ana <- w123 %>% dplyr::filter(forcedchoice == 1 & forcedchoice_w2 == 0)

ctrl_mean_sd_N <- function(y) {
  z <- ana[[y]][ ana$entertainment == 1 ]
  c(mean = mean(z, na.rm = TRUE), sd = sd(z, na.rm = TRUE), N = sum(!is.na(z)))
}


# --- First PC: control-group means/SD (and Ns) by wave
pc_w1_ctrl <- ctrl_mean_sd_N("pcw1")
pc_w2_ctrl <- ctrl_mean_sd_N("pcw2")
pc_w3_ctrl <- ctrl_mean_sd_N("pcw3")

# Formatters
fmt  <- function(x) sprintf("%.2f", x)
fmt0 <- function(x) format(round(x), big.mark = ",")

# Build caption strings

caption_pc <- paste0(
  "Control (Entertainment) means [SD], N — Wave 1: ", fmt(pc_w1_ctrl["mean"]), " [", fmt(pc_w1_ctrl["sd"]), "], N=", fmt0(pc_w1_ctrl["N"]), "; ",
  "Wave 2: ", fmt(pc_w2_ctrl["mean"]), " [", fmt(pc_w2_ctrl["sd"]), "], N=", fmt0(pc_w2_ctrl["N"]), "; ",
  "Wave 3: ", fmt(pc_w3_ctrl["mean"]), " [", fmt(pc_w3_ctrl["sd"]), "], N=", fmt0(pc_w3_ctrl["N"]), ". "
)

cat("\n[Figure 3 — First PC note]\n",     caption_pc,    "\n", sep = "")


###############################################
# Appendix C - Persuasive effects by subgroups (media preference)
###############################################

w123$pid3w1_leaner <- ifelse(w123$party1 == 1, 'Democrat',
                             ifelse(w123$party1 == 2, 'Republican',
                                    ifelse(w123$party1 %in% c(3, 4) & w123$party4 == 1, 'Republican',
                                           ifelse(w123$party1 %in% c(3, 4) & w123$party4 == 2, 'Democrat',
                                                  ifelse(w123$party1 %in% c(3, 4) & w123$party4 == 3, 'Independent', NA)))))

w123$pid3w1_leaner <- factor(w123$pid3w1_leaner, levels = c('Democrat', 'Republican', 'Independent'))



# Compute the cross table
cross_table <- table(w123$pid3w1_leaner, w123$med_choice)
latex_table <- xtable(cross_table, 
                      caption = "Cross Table of Party ID and Media Preference",
                      label = "tab:cross_table")
print(latex_table, type = "latex", include.rownames = TRUE)
print(latex_table, type = "latex", file = "Output/tableA3_cross_table.tex", include.rownames = TRUE)


# First Principal Component Plot
pc_w1_fox <- lm(pcw1 ~ msnbc + fox, data = filter(w123,forcedchoice==1 & forcedchoice_w2==0 & med_pref == "Fox"))
pc_w2_fox <- lm(pcw2 ~ msnbc + fox, data = filter(w123,forcedchoice==1 & forcedchoice_w2==0 & med_pref == "Fox"))
pc_w3_fox <- lm(pcw3 ~ msnbc + fox, data = filter(w123,forcedchoice==1 & forcedchoice_w2==0 & med_pref == "Fox"))
pc_w1_msnbc <- lm(pcw1 ~ msnbc + fox, data = filter(w123,forcedchoice==1 & forcedchoice_w2==0 & med_pref == "MSNBC"))
pc_w2_msnbc <- lm(pcw2 ~ msnbc + fox, data = filter(w123,forcedchoice==1 & forcedchoice_w2==0 & med_pref == "MSNBC"))
pc_w3_msnbc <- lm(pcw3 ~ msnbc + fox, data = filter(w123,forcedchoice==1 & forcedchoice_w2==0 & med_pref == "MSNBC"))
pc_w1_ent <- lm(pcw1 ~ msnbc + fox, data = filter(w123,forcedchoice==1 & forcedchoice_w2==0  & med_pref == "Entertainment"))
pc_w2_ent <- lm(pcw2 ~ msnbc + fox, data = filter(w123,forcedchoice==1 & forcedchoice_w2==0 & med_pref == "Entertainment"))
pc_w3_ent <- lm(pcw3 ~ msnbc + fox, data = filter(w123,forcedchoice==1 & forcedchoice_w2==0 & med_pref == "Entertainment"))

marijuana_ests_fox <- tidy(pc_w1_fox) %>% mutate(Wave = 1,`Media Preference` = "Stated Preference: Fox") %>%
  bind_rows(tidy(pc_w2_fox) %>% mutate(Wave = 2,`Media Preference` = "Stated Preference: Fox")) %>%
  bind_rows(tidy(pc_w3_fox) %>% mutate(Wave = 3,`Media Preference` = "Stated Preference: Fox"))
marijuana_ests_msnbc <- tidy(pc_w1_msnbc) %>% mutate(Wave = 1,`Media Preference` = "Stated Preference: MSNBC") %>%
  bind_rows(tidy(pc_w2_msnbc) %>% mutate(Wave = 2,`Media Preference` = "Stated Preference: MSNBC")) %>%
  bind_rows(tidy(pc_w3_msnbc) %>% mutate(Wave = 3,`Media Preference` = "Stated Preference: MSNBC"))
marijuana_ests_ent <- tidy(pc_w1_ent) %>% mutate(Wave = 1,`Media Preference` = "Stated Preference: Entertainment") %>%
  bind_rows(tidy(pc_w2_ent) %>% mutate(Wave = 2,`Media Preference` = "Stated Preference: Entertainment")) %>%
  bind_rows(tidy(pc_w3_ent) %>% mutate(Wave = 3,`Media Preference` = "Stated Preference: Entertainment"))

marijuana_ests_all <- marijuana_ests_fox %>% 
  bind_rows(marijuana_ests_msnbc) %>%
  bind_rows(marijuana_ests_ent) %>%
  mutate(`Media Preference` = factor(`Media Preference`,levels=c("Stated Preference: MSNBC","Stated Preference: Entertainment","Stated Preference: Fox"),ordered=T),
         upper = estimate + qnorm(0.975)*std.error,
         lower = estimate + qnorm(0.025)*std.error,
         upper_90 = estimate + qnorm(0.95)*std.error,
         lower_90 = estimate + qnorm(0.05)*std.error,
         term = dplyr::recode(term,"msnbc"="MSNBC","fox" = "Fox"))

labs <- tibble(
  Wave = c(1.5, 1.5), 
  estimate = c(0.4, -0.3),
  label = c("Fox rather\nthan Entertainment", "MSNBC rather\nthan Entertainment"),
  term = c("Fox", "MSNBC"),
  `Media Preference` = "Stated Preference: MSNBC"
) %>%
  mutate(`Media Preference` = factor(`Media Preference`,
                                     levels = c("Stated Preference: MSNBC", "Stated Preference: Entertainment", "Stated Preference: Fox"),
                                     ordered = T))

pdf('Output/figA1_pc.pdf', width=6, height=8)
print(ggplot(filter(marijuana_ests_all,term!="(Intercept)"), aes(y=estimate, x=Wave, color=term)) +
  geom_hline(yintercept = 0,lty=2) + 
  geom_point(position=position_dodge(width=0.25))+ 
  geom_line(aes(group=term),position=position_dodge(width=0.25)) + 
  geom_errorbar(aes(ymin=lower, ymax=upper,width=0),size=0.5,position=position_dodge(width=0.25)) + 
  geom_errorbar(aes(ymin=lower_90, ymax=upper_90,width=0),size=1,position=position_dodge(width=0.25)) + 
  facet_wrap(~`Media Preference`,ncol=1) + 
  geom_text(data=labs,aes(label=label)) +
  theme_bw() + 
  scale_x_continuous('Survey Wave',breaks=c(1,2,3),limits=c(0.75,3.25)) + 
  scale_y_continuous('Effect of Treatment on\nFirst Principal Component',
                     breaks=seq(-0.6,0.6,0.3),
                     labels=c("\n\n-0.6\n(More\nliberal)", "-0.3",  "0.0","0.3", "(More\nconservative)\n0.6\n\n"),limits=c(-0.72,0.72)) + 
  scale_color_manual("Treatment",breaks = c("Fox","MSNBC"),values=c(red_mit,blue_mit)) +
  theme(legend.position="none",strip.background = element_rect(fill="lightgrey"),axis.title.y = element_text(margin = margin(r = 0))))
dev.off()



###############################################
# Figure 6 - Accumulation effects
###############################################
w123 <- w123 %>%
  mutate(dosage = case_when(msnbc_w2==1 & msnbc==1 ~ "MSNBC twice",
                            msnbc_w2==1 & entertainment==1 ~ "MSNBC once",
                            fox_w2==1 & fox==1 ~ "Fox twice",
                            fox_w2==1 & entertainment==1 ~ "Fox once",
                            entertainment_w2==1 & entertainment==1 ~ "Entertainment twice"
  ),
  dosage = relevel(factor(dosage),ref = "Entertainment twice"))


# First Principal Component Plot
pc_w2_omnibus <- lm(pcw2 ~ dosage, data = filter(w123,forcedchoice==1 & forcedchoice_w2==1))

marijuana_ests <- tidy(pc_w2_omnibus) %>%
  mutate(upper = estimate + qnorm(0.975)*std.error,
         lower = estimate + qnorm(0.025)*std.error,
         upper_90 = estimate + qnorm(0.95)*std.error,
         lower_90 = estimate + qnorm(0.05)*std.error,
         term = str_replace(term,"dosage",""),
         dosage = str_replace(term,"Fox ",""),
         dosage = str_replace(dosage,"MSNBC ",""),
         term = str_replace(term," once",""),
         term = str_replace(term," twice",""),
         term = factor(term,levels=c("(Intercept)","Fox","MSNBC"),ordered=T)
  )

pdf('Output/fig6_pc.pdf', width=6, height=4)
print(ggplot(filter(marijuana_ests,term != "(Intercept)"), aes(y=estimate, color=term, alpha=dosage,shape=dosage, x=1)) +
  geom_hline(yintercept = 0,lty=2) + 
  geom_point(size=3,position=position_dodge(width=0.2)) + 
  # geom_line(aes(group=term),position=position_dodge(width=0.25)) + 
  geom_errorbar(aes(ymin=lower, ymax=upper,width=0),size=0.5,position=position_dodge(width=0.2)) + 
  geom_errorbar(aes(ymin=lower_90, ymax=upper_90,width=0),size=1,position=position_dodge(width=0.2)) + 
  facet_wrap(~term,scales = "free_x") + 
  # geom_text(data=labs,aes(label=label)) +
  theme_bw() + 
  scale_x_continuous('Waves Treated with Partisan Media',breaks=c(0.95,1.05),labels=c("Wave 2\nonly","Wave 1\nand 2")) + 
  scale_y_continuous('Effect of Treatment on\nFirst Principal Component',
                     breaks=seq(-0.5,0.5,0.25),
                     labels=c("\n\n-0.5\n(More\nliberal)", "-0.25",  "0.0",  "0.25",  "(More\nconservative)\n0.5\n\n"),limits=c(-0.51,0.51)) + 
  scale_color_manual("Treatment",breaks = c("Fox","MSNBC"),values=c(red_mit,blue_mit)) +
  scale_alpha_discrete("Wave 1 Treatment",range=c(0.25,3)) + 
  scale_shape_manual("Wave 1 Treatment",values=c("once"=17,"twice"=16)) + 
  theme(legend.position="none",strip.background = element_rect(fill="lightgrey"),axis.title.y = element_text(margin = margin(r = -20))))
dev.off()


# Effect ratio estimates 

save_ratio_tex <- function(ratio, filename) {
  formatted_ratio <- format(round(ratio, 2))
  writeLines(formatted_ratio, paste0("Output/Estimates/", filename, ".tex"))
}
coef_summary <- summary(pc_w2_omnibus)$coefficients

fox_once_effect <- coef_summary["dosageFox once", "Estimate"]
fox_twice_effect <- coef_summary["dosageFox twice", "Estimate"]
fox_twice_se <- coef_summary["dosageFox twice", "Std. Error"]
fox_ratio_lower <- (fox_twice_effect - 1.96 * fox_twice_se) / fox_once_effect
fox_ratio_upper <- (fox_twice_effect + 1.96 * fox_twice_se) / fox_once_effect
save_ratio_tex(fox_ratio_lower, "fox_two_dose_ratio_lower")
save_ratio_tex(fox_ratio_upper, "fox_two_dose_ratio_upper")


msnbc_once_effect <- coef_summary["dosageMSNBC once", "Estimate"]
msnbc_twice_effect <- coef_summary["dosageMSNBC twice", "Estimate"]
msnbc_twice_se <- coef_summary["dosageMSNBC twice", "Std. Error"]
msnbc_ratio_lower <- (msnbc_twice_effect - 1.96 * msnbc_twice_se) / msnbc_once_effect
msnbc_ratio_upper <- (msnbc_twice_effect + 1.96 * msnbc_twice_se) / msnbc_once_effect
save_ratio_tex(msnbc_ratio_lower, "msnbc_two_dose_ratio_lower")
save_ratio_tex(msnbc_ratio_upper, "msnbc_two_dose_ratio_upper")

print(paste("Fox ratio range:", round(fox_ratio_lower, 3), "to", round(fox_ratio_upper, 3)))
print(paste("MSNBC ratio range:", round(msnbc_ratio_lower, 3), "to", round(msnbc_ratio_upper, 3)))


###############################################
# Figure 7 - Minimum Sample Size
###############################################

calc_mss <- function(effect, power=.8, attrition=0.12, v, alpha=.05, frac_used=1/9){
  M <- qnorm(1-alpha/2, lower.tail=TRUE) + qnorm(power, lower.tail = TRUE)
  (M *  sqrt(2*v/((1-attrition) * .25*(frac_used))))^2/effect^2
}

make_mss_df <- function(baseline_multiple, multiples, frac_used){
  mss <- data.frame(multiples, 
                    Fox = unlist(lapply(effects_fox*baseline_multiple, calc_mss, v=var(w123$pcw1, na.rm=TRUE), frac_used=frac_used)),
                    Fox_lower = unlist(lapply(upper_fox*baseline_multiple, calc_mss, v=var(w123$pcw1, na.rm=TRUE), frac_used=frac_used)),
                    MSNBC = unlist(lapply(effects_msnbc*baseline_multiple, calc_mss, v=var(w123$pcw1, na.rm=TRUE), frac_used=frac_used)),
                    MSNBC_lower = unlist(lapply(upper_msnbc*baseline_multiple, calc_mss, v=var(w123$pcw1, na.rm=TRUE), frac_used=frac_used)))
  
  mss$multiples <- mss$multiples + 1
  mss$baseline_multiple <- baseline_multiple
  mss$frac_used <- frac_used
  mss
}

make_mss_df_base_effect <- function(frac_used, baseline_multiples, multiples){
  do.call(rbind, lapply(baseline_multiples, make_mss_df, multiples= multiples, frac_used=frac_used))
}

mod <- lm(pcw1 ~ fox + msnbc, data = filter(w123, forcedchoice == 1 & forcedchoice_w2 == 0))
multiples <- seq(0, 3, by=.05)

effects_msnbc <- rep(coef(mod)[["msnbc"]], length(multiples)) * multiples  
upper_msnbc <- rep(coef(mod)[["msnbc"]], length(multiples)) * multiples - 1.96 * sqrt(diag(vcov(mod))["msnbc"])

effects_fox <- rep(coef(mod)[["fox"]], length(multiples)) * multiples  
upper_fox <- rep(coef(mod)[["fox"]], length(multiples)) * multiples + 1.96 * sqrt(diag(vcov(mod))["fox"])

fracs <- c(1/9, 1/5, 1/3)
mss <- do.call(rbind, lapply(fracs, make_mss_df_base_effect, baseline_multiples=1:2, multiples=multiples))
mss$n_groups <- paste0(mss$frac_used^(-1), ' Treatment Groups')

names(mss)[names(mss) == "multiples"] <- "Gamma"

mss_long <- pivot_longer(mss[, c('Gamma', 'Fox', "MSNBC", "baseline_multiple", "n_groups")], 
                         cols=c('Fox','MSNBC'), 
                         names_to = "name", 
                         values_to = "value")
mss_long$baseline_multiple <- ifelse(mss_long$baseline_multiple == 1, 'Empirical Effect', '2X Empirical Effect')

pdf('Output/fig7_mss.pdf', width=10, height=5)
print(ggplot(mss_long[mss_long$Gamma <= 2.5 & mss_long$Gamma > 1.50 & mss_long$baseline_multiple=='Empirical Effect',], 
       aes(x=Gamma, y=value, color=n_groups)) + 
  geom_point() + 
  geom_line() + 
  geom_hline(yintercept=4886, alpha=.5) + 
  geom_hline(yintercept=0, alpha=.25, linetype=2) +
  facet_grid(~name) +
  theme_bw() + 
  theme(legend.title = element_blank(),
        strip.background = element_rect(fill="lightgrey"), 
        axis.title.y = element_text(margin = margin(r = 20)),
        plot.background = element_rect(fill = 'white'), 
        panel.background = element_rect(fill = 'white')) + 
  scale_x_continuous(labels = scales::percent) +
  scale_y_continuous(labels = scales::comma) +
  labs(x = 'Effect of Two Doses as % of Effect of One',
       y = 'Minimum Sample Size (MSS)'))
dev.off()


###############################################
# Appendix B - Power Calculation
###############################################

table_data <- mss %>%
  filter(baseline_multiple == 1, 
         n_groups == "9 Treatment Groups", 
         Gamma %in% seq(1.25, 4, by = 0.25))
table_data <- table_data[, c("Gamma", "Fox", "Fox_lower", "MSNBC", "MSNBC_lower")]

table_data <- table_data %>%
  mutate(across(c(Fox, Fox_lower, MSNBC, MSNBC_lower), round))

table_a1 <- kable(table_data, 
                  format = "latex",
                  col.names = c("$\\gamma$", "MSS", "95\\% CI Lower Bound", "MSS", "95\\% CI Lower Bound"),
                  caption = "Minimum Sample Size (MSS) Needed for 80\\% Power By Proportional Change in Treatment Effect ($\\gamma$) Based on Wave 1 Treatment",
                  align = c('c', 'r', 'r', 'r', 'r'),
                  booktabs = TRUE,
                  escape = FALSE) %>%
  kable_styling(latex_options = c("striped", "hold_position")) %>%
  add_header_above(c(" " = 1, "Fox" = 2, "MSNBC" = 2))

writeLines(table_a1, "Output/tableA2.tex")

###############################################
# Save estimates
###############################################

save_estimate_tex <- function(estimate, filename) {
  formatted_estimate <- format(round(estimate, 3))
  writeLines(formatted_estimate, paste0("Output/Estimates/", filename, ".tex"))
}

extract_and_save_estimates <- function(model, wave, outcome) {
  coef_data <- coef(summary(model))
  treatments <- c("fox", "msnbc")
  
  for (treatment in treatments) {
    estimate <- coef_data[treatment, "Estimate"]
    std_error <- coef_data[treatment, "Std. Error"]
    save_estimate_tex(estimate, paste0("wave", wave, "_", outcome, "_", treatment, "_effect"))
    save_estimate_tex(estimate - 1.96 * std_error, paste0("wave", wave, "_", outcome, "_", treatment, "_effect_lower"))
    save_estimate_tex(estimate + 1.96 * std_error, paste0("wave", wave, "_", outcome, "_", treatment, "_effect_upper"))
  }
}

extract_and_save_estimates(pc_w1, 1, "marijuana_pc")
extract_and_save_estimates(pc_w2, 2, "marijuana_pc")
extract_and_save_estimates(pc_w3, 3, "marijuana_pc")